package scu.maqiang.electromagnetics;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.*;

import java.util.Arrays;
import java.util.HashSet;
import java.util.stream.DoubleStream;
import java.util.stream.IntStream;

public class Waveguide_2D {
    public static void main(String[] args) {
        //FEM2D_wg4_Rectangle();
        //FEM2D_wg4_Circle();
        //FEM2_wg1_Rectangle();
        //FEM_wg1_StepDiscontinuous();
        //FEM_wg1_PlateMountedOnGround();
        //FEM2D_wg1_HornMountedOnGround();
        //FEM2D_wg1_HornFreespace();
        //FEM2D_wg2_RectanglePEC();
        //FEM2D_wg2_CircleDiel();
        //FEM2D_wg3_LineSource_RectanglePEC();
        FEM2D_wg4_Circle();
    }


    public static void FEM2D_wg4_Rectangle() {
        Mesh2T3 mesh = new Mesh2T3().square2D(160, 80, T3Type.Left);
        mesh.scale(2.0, 1.0);
        mesh.translate(-1.0, -0.5);
        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SRMatrix B = new SRMatrix(fs.GetNdof());
        fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
        fs.assembleMass(new double[]{1.0}, BVPType.CONSIST_MASS, B);
        double[] x = MVO.rand(fs.GetNdof());
//        fs.applyBC_MBN(A, x, Direct.All, 0.0, 1, 2, 3, 4);

        int nEV = 18;
        double[] eigenvalues = new double[nEV];
        double[][] eigenvectors = new double[nEV][fs.GetNdof()];
        EigenSSolver eigenSolver = new EigenSSolver();
        double alpha = 10.0;
        eigenSolver.inversePowerTran_PCGSSOR(A, B, x, alpha, nEV, eigenvalues, eigenvectors);
        //DoubleStream.of(eigenvalues).forEach(System.out::println);
        mesh.toTecplot("EigenModeRectangleTE.dat", eigenvectors);
        double[] kc = new double[nEV];
        Arrays.setAll(kc, i->Math.sqrt(eigenvalues[i]));
        DoubleStream.of(kc).forEach(System.out::println);
    }

    public static void FEM2D_wg4_Circle() {
        String meshFile = "CircleMesh.mmsh";
        Mesh2T3 mesh = new Mesh2T3().readMesh(meshFile, MeshFormat.mmsh);
        mesh.saveMesh(MeshFormat.msh, "CircleMesh.msh");
        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SRMatrix B = new SRMatrix(fs.GetNdof());
        fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
        fs.assembleMass(new double[]{1.0}, BVPType.CONSIST_MASS, B);
        double[] x = MVO.rand(fs.GetNdof());
        fs.applyBC_MBN(A, x, Direct.All, 0.0, 0);

        int nEV = 10;
        double[] eigenvalues = new double[nEV];
        double[][] eigenvectors = new double[nEV][fs.GetNdof()];
        EigenSSolver eigenSolver = new EigenSSolver();
        double alpha = 10.0;
        //eigenSolver.inversePowerTran_PCGSSOR(A, B, x, alpha, nEV, eigenvalues, eigenvectors);
        eigenSolver.inversePowerTran_DSS(A, B, x, alpha, nEV, eigenvalues, eigenvectors);
        //DoubleStream.of(eigenvalues).forEach(System.out::println);
        mesh.toTecplot("EigenModeCircleTE.dat", eigenvectors);
        double[] kc = new double[nEV];
        Arrays.setAll(kc, i->Math.sqrt(eigenvalues[i]));
        DoubleStream.of(kc).forEach(System.out::println);
    }

    public static void FEM2_wg1_Rectangle() {

        //参数设置
        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 2000;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency

        Mesh2T3 mesh = new Mesh2T3().square2D(220, 40, T3Type.Left);
        mesh.scale(0.55, 0.1);
        int pmlDomainLabel = 10;
        ScalarFunc pmlDomainFunc = (xy, label, param) -> xy[0] > 0.5?1.0:-1.0;
        mesh.setDomainLabel(pmlDomainFunc, null, pmlDomainLabel);

        int interfaceBoundaryLabel  = 20;
        ScalarFunc pmlInnerBoundary = (xy, label, param) -> Math.abs(xy[0] - 0.5) < 1e-5?1.0:-1.0;
        mesh.addInterfaceBoundary(pmlInnerBoundary, null, interfaceBoundaryLabel);

        int pmlOuterBoundaryLabel = 2;

        double dx = 0.1;
        int nmode = 1;
        //ScalarCFunc boundaryFunc = (xy, label, param) -> new Complex(Math.sin(nmode*Math.PI*xy[1]/dx), 0);
        ScalarCFunc boundaryFunc = (xy, label, param) -> new Complex(Math.cos(nmode*Math.PI*xy[1]/dx), 0);
        FES2T31 fs = new FES2T31(mesh);
        SCMatrix A = new SCMatrix(fs.GetNdof());

        double mur = 1.0;
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, 0);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, 0);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundaryLabel, pmlOuterBoundaryLabel, pmlDomainLabel);

        //A.sort();
        //System.out.println(A.display());

        Complex[] RHS = Complex.zeros(fs.GetNdof());
        //fs.applyBC_MBN(A, Direct.All, 1, 2, 3);
        //fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), 1, 2, 3);

        fs.applyBC_MBN(A, Direct.All, 4);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{boundaryFunc}, 4);

        IterSCSolver solver = new IterSCSolver(A);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 2);


        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[3][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        mesh.toTecplot("Waveguide_rectangle.dat", result);
    }

    public static void FEM_wg1_StepDiscontinuous() {

        //参数设置
        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 4000;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency

        String meshFile = "StepDiscontinuousMesh.mmsh";
        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(meshFile, MeshFormat.mmsh);
        int pmlDomainLabel = 10;
        ScalarFunc pmlDomainFunc = (xy, label, param) -> xy[0] > 0.5?1.0:-1.0;
        mesh.setDomainLabel(pmlDomainFunc, null, pmlDomainLabel);

        int interfaceBoundaryLabel  = 20;
        ScalarFunc pmlInnerBoundary = (xy, label, param) -> Math.abs(xy[0] - 0.5) < 1e-5?1.0:-1.0;
        mesh.addInterfaceBoundary(pmlInnerBoundary, null, interfaceBoundaryLabel);

        int pmlOuterBoundaryLabel = 2;
        ScalarFunc pmlOuterBoundary = (xy, label, param) -> Math.abs(xy[0] - 0.55) < 1e-4?1.0:-1.0;
        mesh.setBoundaryLabel(pmlOuterBoundary, null, pmlOuterBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(pmlOuterBoundaryLabel).size());

        int leftBoundaryLabel = 4;
        ScalarFunc leftBoundary = (xy, label, param) -> Math.abs(xy[0]) < 1e-5?1.0:-1.0;
        mesh.setBoundaryLabel(leftBoundary, null, leftBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(leftBoundaryLabel).size());

        double dx = 0.1;
        int nmode = 1;
        //ScalarCFunc boundaryFunc = (xy, label, param) -> new Complex(Math.sin(nmode*Math.PI*xy[1]/dx), 0);
        ScalarCFunc boundaryFunc = (xy, label, param) -> new Complex(Math.cos(nmode*Math.PI*xy[1]/dx), 0);
        FES2T31 fs = new FES2T31(mesh);
        SCMatrix A = new SCMatrix(fs.GetNdof());

        double mur = 1.0;
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, 0);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, 0);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundaryLabel, pmlOuterBoundaryLabel, pmlDomainLabel);

        //A.sort();
        //System.out.println(A.display());

        Complex[] RHS = Complex.zeros(fs.GetNdof());
        //fs.applyBC_MBN(A, Direct.All, 0, 2);
        //fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), 0, 2);

        fs.applyBC_MBN(A, Direct.All, 4);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{boundaryFunc}, 4);

        IterSCSolver solver = new IterSCSolver(A);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 2);


        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[3][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        mesh.toTecplot("Waveguide_StepDiscontinuous.dat", result);
    }

    public static void FEM_wg1_PlateMountedOnGround() {
        //参数设置
        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 2000;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency

        String meshFile = "PlateMountedOnGround.mmsh";
        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(meshFile, MeshFormat.mmsh);
        int pmlDomainLabel = 0;
        int freeSpaceDomainLabel = 10;
        ScalarFunc freeSpaceDomainFunc = (xy, label, param) -> Math.abs(xy[0] - 0.1) < 0.1 && Math.abs(xy[1] - 0.05) < 0.15?1.0:-1.0;
        mesh.setDomainLabel(freeSpaceDomainFunc, null, freeSpaceDomainLabel);

        int interfaceBoundaryLabel  = 20;
        ScalarFunc pmlInnerBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[1] - 0.05)-0.15) < 1.0e-4 && Math.abs(xy[0] - 0.15) < 0.05) {
                return 1.0;
            }
            if(Math.abs(Math.abs(xy[0] - 0.2)) < 1.0e-4 && Math.abs(xy[1] - 0.05) < 0.15) {
                return 1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(pmlInnerBoundary, null, interfaceBoundaryLabel);

        int pmlOuterBoundaryLabel = 2;
        ScalarFunc pmlOuterBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[1] - 0.05)-0.2) < 1.0e-4 && Math.abs(xy[0] - 0.175) < 0.75) {
                return 1.0;
            }
            if(Math.abs(Math.abs(xy[0] - 0.25)) < 1.0e-4) {
                return 1.0;
            }
            return -1.0;
        };
        mesh.setBoundaryLabel(pmlOuterBoundary, null, pmlOuterBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(pmlOuterBoundaryLabel).size());

        int leftBoundaryLabel = 4;
        ScalarFunc leftBoundary = (xy, label, param) -> Math.abs(xy[0]) < 1e-5?1.0:-1.0;
        mesh.setBoundaryLabel(leftBoundary, null, leftBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(leftBoundaryLabel).size());

        double dx = 0.1;
        int nmode = 1;
        ScalarCFunc boundaryFunc = (xy, label, param) -> new Complex(Math.sin(nmode*Math.PI*xy[1]/dx), 0);
        //ScalarCFunc boundaryFunc = (xy, label, param) -> new Complex(Math.cos(nmode*Math.PI*xy[1]/dx), 0);
        FES2T31 fs = new FES2T31(mesh);
        SCMatrix A = new SCMatrix(fs.GetNdof());

        double mur = 1.0;
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, freeSpaceDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, freeSpaceDomainLabel);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundaryLabel, pmlOuterBoundaryLabel, pmlDomainLabel);

        //A.sort();
        //System.out.println(A.display());

        Complex[] RHS = Complex.zeros(fs.GetNdof());
        fs.applyBC_MBN(A, Direct.All, 0);
        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), 0);

        fs.applyBC_MBN(A, Direct.All, 4);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{boundaryFunc}, 4);

        IterSCSolver solver = new IterSCSolver(A);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 2);


        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[3][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        mesh.toTecplot("Waveguide_PlateMountedOnGround.dat", result);

    }

    public static void FEM2D_wg1_HornMountedOnGround() {
        //参数设置
        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 2000;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency

        String meshFile = "HornMountedOnGround.mmsh";
        Mesh2T3 mesh = new Mesh2T3().readMesh(meshFile, MeshFormat.mmsh);
        int pmlDomainLabel = 0;
        int freeSpaceDomainLabel = 10;
        ScalarFunc freeSpaceDomainFunc = (xy, label, param) -> Math.abs(xy[0] - 0.125) < 0.125 && Math.abs(xy[1] - 0.05) < 0.2?1.0:-1.0;
        mesh.setDomainLabel(freeSpaceDomainFunc, null, freeSpaceDomainLabel);

        int interfaceBoundaryLabel  = 20;
        ScalarFunc pmlInnerBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[1] - 0.05) - 0.2) < 1.0e-4 && Math.abs(xy[0] - 0.2) < 0.05) {
                return 1.0;
            }
            if(Math.abs(Math.abs(xy[0] - 0.25)) < 1.0e-4 && Math.abs(xy[1] - 0.05) < 0.2) {
                return 1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(pmlInnerBoundary, null, interfaceBoundaryLabel);

        int pmlOuterBoundaryLabel = 2;
        ScalarFunc pmlOuterBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[1] - 0.05)-0.25) < 1.0e-4 && Math.abs(xy[0] - 0.225) < 0.75) {
                return 1.0;
            }
            if(Math.abs(Math.abs(xy[0] - 0.3)) < 1.0e-4) {
                return 1.0;
            }
            return -1.0;
        };
        mesh.setBoundaryLabel(pmlOuterBoundary, null, pmlOuterBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(pmlOuterBoundaryLabel).size());

        int leftBoundaryLabel = 4;
        ScalarFunc leftBoundary = (xy, label, param) -> Math.abs(xy[0]) < 1e-5?1.0:-1.0;
        mesh.setBoundaryLabel(leftBoundary, null, leftBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(leftBoundaryLabel).size());

        double dx = 0.1;
        int nmode = 1;
        ScalarCFunc boundaryFunc = (xy, label, param) -> new Complex(Math.sin(nmode*Math.PI*xy[1]/dx), 0);
        //ScalarCFunc boundaryFunc = (xy, label, param) -> new Complex(Math.cos(nmode*Math.PI*xy[1]/dx), 0);
        FES2T31 fs = new FES2T31(mesh);
        SCMatrix A = new SCMatrix(fs.GetNdof());

        double mur = 1.0;
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, freeSpaceDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, freeSpaceDomainLabel);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundaryLabel, pmlOuterBoundaryLabel, pmlDomainLabel);

        //A.sort();
        //System.out.println(A.display());

        Complex[] RHS = Complex.zeros(fs.GetNdof());
        fs.applyBC_MBN(A, Direct.All, 0);
        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), 0);

        fs.applyBC_MBN(A, Direct.All, 4);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{boundaryFunc}, 4);

        IterSCSolver solver = new IterSCSolver(A);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 2);


        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[3][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        mesh.toTecplot("Waveguide_HornMountedOnGround.dat", result);
    }

    public static void FEM2D_wg1_HornFreespace() {

        //参数设置
        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 2000;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency
        double dx = 0.1;
        int nmode = 1;
        double kc = nmode * Math.PI / dx;
        double fc = kc/(2*Math.PI*Math.sqrt(mu0*e0));
        double beta = Math.sqrt(k0 * k0 - kc * kc);

        String meshFile = "HornFreespace.mmsh";
        Mesh2T3 mesh = new Mesh2T3().readMesh(meshFile, MeshFormat.mmsh);
        int pmlDomainLabel = 0;
        int freeSpaceDomainLabel = 10;
        ScalarFunc freeSpaceDomainFunc = (xy, label, param) -> Math.abs(xy[0] - 0.175) < 0.125 && Math.abs(xy[1] - 0.05) < 0.2?1.0:-1.0;
        mesh.setDomainLabel(freeSpaceDomainFunc, null, freeSpaceDomainLabel);

        int interfaceBoundaryLabel  = 20;
        ScalarFunc pmlInnerBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[1] - 0.05) - 0.2) < 1.0e-4 && Math.abs(xy[0] - 0.175) < 0.125) {
                return 1.0;
            }
            if(Math.abs(Math.abs(xy[0] - 0.175) - 0.125) < 1.0e-4 && Math.abs(xy[1] - 0.05) < 0.2) {
                return 1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(pmlInnerBoundary, null, interfaceBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(interfaceBoundaryLabel).size());

        int pmlOuterBoundaryLabel = 0;
        System.out.println(mesh.extractBoundaryNodes(pmlOuterBoundaryLabel).size());

        int leftBoundaryLabel = 4;
        ScalarFunc leftBoundary = (xy, label, param) -> Math.abs(xy[0] - 0.055) < 1e-4 && Math.abs(xy[1] - 0.05) < 0.05?1.0:-1.0;
        mesh.addInterfaceBoundary(leftBoundary, null, leftBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(leftBoundaryLabel).size());

        int zeroBoundaryLabel = 30;
        ScalarFunc zeroBoundary = (xy, label, param) -> {
            if(Math.abs(xy[0] - 0.05) < 0.05 && Math.abs(Math.abs(xy[1] - 0.05) - 0.05) < 1.0e-5) {
                return 1.0;
            }
            if(Math.abs(xy[0] - 0.15) < 0.05) {
                double temp1 = 0.5 * (xy[0] - 0.1) + 0.1;
                double temp2 = -0.5 * (xy[0] - 0.1);
                if(Math.abs(temp1 - xy[1]) < 1.0e-4 || Math.abs(temp2 - xy[1]) < 1.0e-4) {
                    return 1.0;
                }
            }
            return -1.0;
        };

        mesh.addInterfaceBoundary(zeroBoundary, null, zeroBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(zeroBoundaryLabel).size());


        ScalarCFunc boundaryFunc = (xy, label, param) -> {
            double temp = Math.sin(nmode*Math.PI*xy[1]/dx);
            return new Complex(temp * Math.cos(beta * xy[0]), -temp * Math.sin(beta * xy[1]));
        };
        //ScalarCFunc boundaryFunc = (xy, label, param) -> new Complex(Math.cos(nmode*Math.PI*xy[1]/dx), 0);
        FES2T31 fs = new FES2T31(mesh);
        SCMatrix A = new SCMatrix(fs.GetNdof());

        double mur = 1.0;
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, freeSpaceDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, freeSpaceDomainLabel);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundaryLabel, pmlOuterBoundaryLabel, pmlDomainLabel);

        //A.sort();
        //System.out.println(A.display());

        Complex[] RHS = Complex.zeros(fs.GetNdof());
        fs.applyBC_MBN(A, Direct.All, zeroBoundaryLabel);
        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), zeroBoundaryLabel);

        fs.applyBC_MBN(A, Direct.All, leftBoundaryLabel);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{boundaryFunc}, leftBoundaryLabel);

        IterSCSolver solver = new IterSCSolver(A);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 1);


        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[3][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        mesh.toTecplot("Waveguide_HornFreeSpace.dat", result);
    }

    public static void FEM2D_wg2_RectanglePEC() {

        double c0  = 3*1e8;           // m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;     // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        //Input parameters
        //*****************
        double freq    = 2000;         // MHz, frequency
        freq    = freq*1e6;            // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0      = 2 * Math.PI / lambda0; // 1/meter, wavenumber
        double omg     = 2 * Math.PI * freq;           // rad/sec, radial frequency


        double er = 4;  // dielectric constant of dielectric object (used if scat_type='diel')
        double mur = 1; // relative permittivity of object (used if scat_type='diel')

        String meshFile = "Waveguide_RectanglePEC.mmsh";
        Mesh2T3 mesh = new Mesh2T3().readMesh(meshFile, MeshFormat.mmsh);

        //确定自由空间区域与PML区域
        int freeSpaceDomainLabel = 0;

        int pmlDomainLabel_1 = 11;
        ScalarFunc pmlDomain1 = (xy, label, param) -> xy[0] < 0.0? 1.0: -1.0;
        mesh.setDomainLabel(pmlDomain1, null, pmlDomainLabel_1);
        int pmlDomainLabel_2 = 12;
        ScalarFunc pmlDomain2 = (xy, label, param) -> xy[0] > 0.5?1.0:-1.0;
        mesh.setDomainLabel(pmlDomain2, null, pmlDomainLabel_2);

        //设定边界条件
        int waveguideWallLabel = 0;
        int leftBoundary = 4;
        ScalarFunc leftBoundaryFunc = (xy, label, param) -> xy[0] < -0.05?1.0:-1.0;
        mesh.setBoundaryLabel(leftBoundaryFunc, null, leftBoundary);
        System.out.println(mesh.extractBoundaryNodes(leftBoundary).size());
        int rightBoundary = 2;
        ScalarFunc rightBoundaryFunc = (xy, label, param) -> xy[0] > 0.55?1.0:-1.0;
        mesh.setBoundaryLabel(rightBoundaryFunc, null, rightBoundary);
        System.out.println(mesh.extractBoundaryNodes(rightBoundary).size());

        int pecBoundaryLabel = 8;
        ScalarFunc pecBoundaryFunc = (xy, label, param) -> Math.abs(xy[1] - 0.05) < 0.03 && Math.abs(xy[0] - 0.25) < 0.05?1.0:-1.0;
        mesh.setBoundaryLabel(pecBoundaryFunc, null, pecBoundaryLabel);

        int interfaceBoundary_1 = 44;
        ScalarFunc interfactFunc_1 = (xy, label, param) -> Math.abs(xy[0]) < 1.0e-5?1.0:-1.0;
        mesh.addInterfaceBoundary(interfactFunc_1, null, interfaceBoundary_1);

        int interfaceBoundary_2 = 22;
        ScalarFunc interfactFunc_2 = (xy, label, param) -> Math.abs(xy[0] - 0.5) < 1.0e-5?1.0:-1.0;
        mesh.addInterfaceBoundary(interfactFunc_2, null, interfaceBoundary_2);

        double nmode = 1;
        double dx = 0.1;
        double kc = nmode * Math.PI / dx;
        double fc = kc/(2*Math.PI*Math.sqrt(mu0*e0));
        double beta = Math.sqrt(k0 * k0 - kc * kc);

        ScalarCFunc boundaryFunc = (xy, label, param) -> {
            double temp = Math.cos(nmode*Math.PI*xy[1]/dx);
            return new Complex(-temp * Math.cos(beta * xy[0]), +temp * Math.sin(beta * xy[0]));
        };


        FES2T31 fs = new FES2T31(mesh);
        SCMatrix A = new SCMatrix(fs.GetNdof());

        //在自由空间组装整体矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, freeSpaceDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, freeSpaceDomainLabel);

        //在PML区域组装PML矩阵
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundary_1, leftBoundary, pmlDomainLabel_1);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundary_2, rightBoundary, pmlDomainLabel_2);
        //整体矩阵
//        A.sort();
//        System.out.println(A.display());

        //右端项
        Complex[] RHS = Complex.zeros(fs.GetNdof());

        //整个外边界设为0值
        fs.applyBC_MBN(A, Direct.All, waveguideWallLabel, leftBoundary, rightBoundary);
        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), waveguideWallLabel, leftBoundary, rightBoundary);

        //PEC边界，结构区域内部边界
        fs.applyBC_MBN(A, Direct.All, pecBoundaryLabel);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{boundaryFunc}, pecBoundaryLabel);


        IterSCSolver solver = new IterSCSolver(A);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 1);

        //入射场负值
        Complex[] uint = mesh.functionInMesh(boundaryFunc, null);

        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[6][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        Complex[] utot = Complex.add(x, -1.0, uint);
        result[3] = Complex.Real(uint);
        result[4] = Complex.Image(uint);
        result[5] = Complex.abs(uint);
        mesh.toTecplot("Waveguild_RectanglePEC.dat", result);

    }

    public static void FEM2D_wg2_CircleDiel() {

        double c0  = 3*1e8;           // m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;     // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        //Input parameters
        //*****************
        double freq    = 2000;         // MHz, frequency
        freq    = freq*1e6;            // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0      = 2 * Math.PI / lambda0; // 1/meter, wavenumber
        double omg     = 2 * Math.PI * freq;           // rad/sec, radial frequency


        double er = 4;  // dielectric constant of dielectric object (used if scat_type='diel')
        double mur = 1; // relative permittivity of object (used if scat_type='diel')

        String meshFile = "WaveguideR_CircleDiel.mmsh";
        Mesh2T3 mesh = new Mesh2T3().readMesh(meshFile, MeshFormat.mmsh);

        //确定自由空间区域与PML区域
        int freeSpaceDomainLabel = 0;

        int pmlDomainLabel_1 = 11;
        ScalarFunc pmlDomain1 = (xy, label, param) -> xy[0] < 0.0? 1.0: -1.0;
        mesh.setDomainLabel(pmlDomain1, null, pmlDomainLabel_1);
        int pmlDomainLabel_2 = 12;
        ScalarFunc pmlDomain2 = (xy, label, param) -> xy[0] > 0.5?1.0:-1.0;
        mesh.setDomainLabel(pmlDomain2, null, pmlDomainLabel_2);

        String scattererElement = "Waveguide_CircleDiel.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int dielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, dielDomainLabel));
        System.out.println(mesh.extractDomainNodes(dielDomainLabel).size());

        //设定边界条件
        int waveguideWallLabel = 0;
        int leftBoundary = 4;
        ScalarFunc leftBoundaryFunc = (xy, label, param) -> xy[0] < -0.05?1.0:-1.0;
        mesh.setBoundaryLabel(leftBoundaryFunc, null, leftBoundary);
        System.out.println(mesh.extractBoundaryNodes(leftBoundary).size());
        int rightBoundary = 2;
        ScalarFunc rightBoundaryFunc = (xy, label, param) -> xy[0] > 0.55?1.0:-1.0;
        mesh.setBoundaryLabel(rightBoundaryFunc, null, rightBoundary);
        System.out.println(mesh.extractBoundaryNodes(rightBoundary).size());

        int interfaceBoundary_1 = 44;
        ScalarFunc interfactFunc_1 = (xy, label, param) -> Math.abs(xy[0]) < 1.0e-5?1.0:-1.0;
        mesh.addInterfaceBoundary(interfactFunc_1, null, interfaceBoundary_1);

        int interfaceBoundary_2 = 22;
        ScalarFunc interfactFunc_2 = (xy, label, param) -> Math.abs(xy[0] - 0.5) < 1.0e-5?1.0:-1.0;
        mesh.addInterfaceBoundary(interfactFunc_2, null, interfaceBoundary_2);

        double nmode = 1;
        double dx = 0.1;
        double kc = nmode * Math.PI / dx;
        double fc = kc/(2*Math.PI*Math.sqrt(mu0*e0));
        double beta = Math.sqrt(k0 * k0 - kc * kc);

        ScalarCFunc boundaryFunc = (xy, label, param) -> {
            double temp = Math.cos(nmode*Math.PI*xy[1]/dx);
            return new Complex(-temp * Math.cos(beta * xy[0]), +temp * Math.sin(beta * xy[0]));
        };

        Complex[] uint = mesh.functionInMesh(boundaryFunc, null);
//        HashSet<Integer> scatternodeSet = mesh.extractDomainNodes(dielDomainLabel);
//        for(int idx: scatternodeSet) {
//            System.out.println(idx + "\t" + uint[idx]);
//        }

        FES2T31 fs = new FES2T31(mesh);
        SCMatrix A = new SCMatrix(fs.GetNdof());

        //在自由空间组装整体矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, freeSpaceDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, freeSpaceDomainLabel);

        //在介电区域组装整体矩阵
        fs.assembleStiff(new double[]{1.0/(mur * 4)}, BVPType.COMMON, A, dielDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0}, BVPType.CONSIST_MASS, A, dielDomainLabel);

        Complex[] RHS = Complex.zeros(fs.GetNdof());
        A.mul(uint, RHS);
        Complex[] newRHS = Complex.zeros(fs.GetNdof());
        HashSet<Integer> nodeSet = mesh.extractDomainNodes(dielDomainLabel);
        for(int idx:nodeSet) {
            newRHS[idx] = RHS[idx];
        }
        //System.out.println(nodeSet.size());
        RHS = newRHS;

        //在PML区域组装PML矩阵
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundary_1, leftBoundary, pmlDomainLabel_1);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundary_2, rightBoundary, pmlDomainLabel_2);

        //Complex.toString(uint);
        //System.out.println(Complex.toString(RHS));
//        A.sort();
//        System.out.println(A.display());
        //System.exit(1);
        //整体矩阵
//        A.sort();
//        System.out.println(A.display());

        //右端项

        //整个外边界设为0值
//        fs.applyBC_MBN(A, Direct.All, waveguideWallLabel);
//        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), waveguideWallLabel);

        //PEC边界，结构区域内部边界
//        fs.applyBC_MBN(A, Direct.All, pecBoundaryLabel);
//        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{boundaryFunc}, pecBoundaryLabel);


        IterSCSolver solver = new IterSCSolver(A);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 1);

        //入射场负值
        //Complex[] uint = mesh.functionInMesh(boundaryFunc, null);

        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[6][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        Complex[] utot = Complex.add(x, -1.0, uint);
        result[3] = Complex.Real(utot);
        result[4] = Complex.Image(utot);
        result[5] = Complex.abs(utot);
        mesh.toTecplot("Waveguide_CircleDiel.dat", result);

    }

    public static void FEM2D_wg3_LineSource_RectanglePEC() {

        //**********
        double c0  = 3*1e8;           // m/sec, velocity of light in free space
        double nu0 = 120 * Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        //Input parameters
        //*****************
        double freq    = 5000;         // MHz, frequency
        freq    = freq*1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0      = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg     = 2*Math.PI*freq;    // rad/sec, radial frequency

        double er = 4;  // dielectric constant of dielectric object (used if scat_type='diel')
        double mur = 1; // relative permittivity of object (used if scat_type='diel')

        double Is = 1e-4;

        String meshFile = "SourceRectangle.mmsh";
        Mesh2T3 mesh = new Mesh2T3().readMesh(meshFile, MeshFormat.mmsh);

        String sourceNodeFile = "SourceNodeInfo.txt";
        String sourceValueFile = "SourceValueInfo.txt";
        int[] sourceNode = MVO.readIntArray(sourceNodeFile);
        Arrays.setAll(sourceNode, i -> sourceNode[i] - 1);
        Complex[] sourceValue = Complex.readArray(sourceValueFile);
//        System.out.println(Complex.toString(sourceValue));
//        System.exit(1);
        //确定自由空间区域与PML区域
        int freeSpaceDomainLabel = 0;

        int pmlDomainLabel_1 = 11;
        ScalarFunc pmlDomain1 = (xy, label, param) -> xy[0] < 0.0? 1.0: -1.0;
        mesh.setDomainLabel(pmlDomain1, null, pmlDomainLabel_1);
        int pmlDomainLabel_2 = 12;
        ScalarFunc pmlDomain2 = (xy, label, param) -> xy[0] > 0.5?1.0:-1.0;
        mesh.setDomainLabel(pmlDomain2, null, pmlDomainLabel_2);

        //设定边界条件
        int waveguideWallLabel = 0;
        int leftBoundary = 4;
        ScalarFunc leftBoundaryFunc = (xy, label, param) -> Math.abs(xy[0] + 0.05) < 1.0e-5?1.0:-1.0;
        mesh.setBoundaryLabel(leftBoundaryFunc, null, leftBoundary);
        System.out.println(mesh.extractBoundaryNodes(leftBoundary).size());
        int rightBoundary = 2;
        ScalarFunc rightBoundaryFunc = (xy, label, param) -> Math.abs(xy[0] - 0.55) < 1.0e-5?1.0:-1.0;
        mesh.setBoundaryLabel(rightBoundaryFunc, null, rightBoundary);
        System.out.println(mesh.extractBoundaryNodes(rightBoundary).size());

        int interfaceBoundary_1 = 44;
        ScalarFunc interfaceFunc_1 = (xy, label, param) -> Math.abs(xy[0]) < 1.0e-5?1.0:-1.0;
        mesh.addInterfaceBoundary(interfaceFunc_1, null, interfaceBoundary_1);

        int interfaceBoundary_2 = 22;
        ScalarFunc interfaceFunc_2 = (xy, label, param) -> Math.abs(xy[0] - 0.5) < 1.0e-5?1.0:-1.0;
        mesh.addInterfaceBoundary(interfaceFunc_2, null, interfaceBoundary_2);

        FES2T31 fs = new FES2T31(mesh);
        SCMatrix A = new SCMatrix(fs.GetNdof());
        //在自由空间组装整体矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, freeSpaceDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, freeSpaceDomainLabel);

        //在PML区域组装PML矩阵
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundary_1, leftBoundary, pmlDomainLabel_1);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A, interfaceBoundary_2, rightBoundary, pmlDomainLabel_2);

        //右端项
        Complex[] RHS = Complex.zeros(fs.GetNdof());

        //整个外边界设为0值
//        fs.applyBC_MBN(A, Direct.All, waveguideWallLabel, leftBoundary, rightBoundary);
//        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), waveguideWallLabel, leftBoundary, rightBoundary);

        //PEC边界，结构区域内部边界
        fs.applyBC_MBN(A, RHS, sourceNode, sourceValue);

        IterSCSolver solver = new IterSCSolver(A);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 2);


        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[3][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        mesh.toTecplot("Waveguild_LineSource_RectanglePEC.dat", result);
    }
}
